function eq = find_agg_mit(eq, param, glob, options)

%finds C and D implied by policies and distributions

tosolve = @(C) kimball_agg(C, eq, param, glob) - glob.agg;
try
   C_out = bisect(tosolve, 1e-12, 1000);
catch
   C_out = 1000000;
   if tosolve(.0001) < -.99
       if tosolve(100000) < -.99
           C_out = 0;
       end
   end
end

q = eq.yf/C_out;
D_out = (sum(eq.L.*upsilonp(q,param).*q))^(-1);

eq.C = C_out;
eq.D = D_out;

end